########### MANUSCRIPT: The Role of Education and Attitudes in Cooking Fuel Choice: Evidence from two states in India
########### JOURNAL: ENERGY FOR SUSTAINABLE DEVELOPMENT
########### AUTHORS: CARLOS F. GOULD/1 AND JOHANNES URPELAINEN/2
########### AFFILIATIONS: 1/COLUMBIA UNIVERSITY MAILMAN SCHOOL OF PUBLIC HEALTH AND 2/JOHNS HOPKINS SCHOOL OF ADVANCED INTERNATIONAL STUDIES
########### PURPOSE: THIS IS CODE #4 THAT MAKES DESCRIPTIVE MAPS

# brief processing goes first --------------------------

# all LPG cylinder costs >1500 rupees per cylinder are considered implausibly large
merged_regs$q127a_1 <- ifelse(merged_regs$q127a_1>1500, NA, merged_regs$q127a_1)

# summarize across districts
districts_costs <- merged_regs %>%
  group_by(q1) %>%
  summarize(lpg_cost = mean(q127a_1, na.rm=T))

# read in survey district data (names, etc.)
districts <- read.csv(file="~/Dropbox/Kerala Rajasthan/Data/Data/In use/Districts.csv")

# merge with the summarize LPG cylinder costs
districts_costs <- districts_costs %>% rename("District_Number" = "q1")
districts <- districts %>%
  left_join(districts_costs)

districts["% lpg"] <- districts$LPG / districts$N
districts$DISTRICT <- as.character(districts$DISTRICT)
districts$ST_NAME <- as.character(districts$ST_NAME)
districts <- districts %>% rename("LPG Cost (INR)" = "lpg_cost")

districts$`LPG Cost (INR)` <- round(districts$`LPG Cost (INR)`, digits = 0)
districts$`% lpg` <- paste0(round(districts$`% lpg`, digits = 2)*100, "%")

# full india shape files
india <- rgdal::readOGR(dsn = "~/Dropbox/LPG Use Survey/Dataverse_Files/", layer = "INDIA")

# data set with districts AND state names (shape files)
districts_full <- rgdal::readOGR(dsn = "~/Dropbox/LPG Use Survey/Dataverse_Files/", layer = "2011_Dist")
districts_full$DISTRICT <- as.character(districts_full$DISTRICT)
india$ST_NAME <- as.character(india$ST_NAME)
districts_full@data <- left_join(districts_full@data, districts, by = "DISTRICT", all = TRUE)
india@data <- left_join(india@data, districts, by = "ST_NAME", all = TRUE)

#remove NA districts
districts_nona <- districts_full[!is.na(districts_full@data$`% lpg`), ]


# FIGURE A3 -------------------------------------------------------------
# Map of study districts showing the fraction of households owning LPG.
# formatting is a bit off from original, but numbers are the same

tm_shape(na.omit(districts_full)) +
  tm_borders(col = "gray50", lwd = 0.5) +
  tm_polygons("% lpg", palette = "Blues") +
  tm_legend(scale = 2.5, aes.color = c(borders = "black"),
            legend.position = c("right", "top"))


# FIGURE A7 -------------------------------------------------------------
# Map showing the geographic distribution of the average price (in rupees) paid by
# households for their most recent LPG cylinder purchase at the district-level. The percent of study
# participants in that district with LPG is overlaid. Prices above 1,500 rupees were removed due to
# implausibility.



tm_shape(na.omit(districts_full)) + 
  tm_polygons("LPG Cost (INR)", palette = "Blues",
              border.alpha = 0,
              showNA = F, colorNA= NULL) +
  tm_shape(na.omit(districts_full)) + 
  tm_text("% lpg", just = "top", 
          showNA = F, textNA= NULL, size = 0.5,
          auto.placement = T) + 
  tm_legend(scale = 2.5, aes.color = c(borders = "black"), 
            legend.position = c("right", "top"))




